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Abstract. Forecasting failure events is one of the most important problems in fracture mechanics and 
related sciences. In this paper, we use the Molchan scheme to investigate the error diagrams in a fracture 
model which has the notable advantage of displaying two completely different regimes according to the 
heterogeneity of the system. In one regime, a characteristic event is observed while for the second regime 
a power-law spectrum of avalanches is obtained reminiscent of self-organized criticality. We find that both 
regimes are different when predicting large avalanches and that, in the second regime, there are non- 
trivial temporal correlations associated to clustering of large events. Finally, we extend the discussion to 
seismology, where both kinds of avalanche size distributions can be seen. 



PACS. 46.50.+a Fracture mechanics, brittleness, fracture and cracks 
62.20.Mk Fatigue, brittleness, fracture, and cracks 



91.45.Vz Fracture and faults 



1 Introduction 

The fracture of heterogeneous materials has been the sub- 
ject of intensive research since many years and is one of the 
oldest concerns of science £Q . This is not only due to the 
evident practical and potential profits, but also because 
the understanding of fracture processes at a basic level 
has shed light on other phenomena like earthquake occur- 
rence • Relatively simple models such as random resis- 
tor networks beam networks pjj, spring-block models 
[uTHTIS] and the so-called Fibre Bundle Models (FBMs) [2 
3,9, rtJlllllll2j have guided our way to more complex mod- 
els of fracture. To date, we actually know many of the 
failure properties and its dynamics. Yet, the range of phe- 
nomena associated with fracture has not been cast into 
a definite physical and theoretical treatment. Some years 
ago, there was a burst of activity in what today we know as 
self-organized criticality 13,14 , a theoretical framework 
that was soon applied to the study of fracture processes 
and avalanche like phenomena in disordered system. Soon 
afterwards, some interesting claimings were put out as, 
for example, the suggestion that the earth's crust is in a 
self-organized critical state, being the Gutenberg-Richter 
(GR) law JH] for the frequency of earthquakes with a 
given magnitude a result of this self-organization process 

EEEIIIB. 

From a more practical point of view, perhaps the most 
important concern related to fracture processes that re- 
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mains unsolved is their predictability 3 , iSj . One is in- 
terested in knowing not only under what circumstances a 
material will fail or how the distribution of energy releases 
looks like, but also when a material will break down or an 
event of a given magnitude (usually we are interested in 
largest quakes) will take place. This problem is far from 
being solved and different approaches have been adopted 
in order to increase the predictive power of actual meth- 
ods | 19) . For instance, in seismicity, there is an unavoidable 
degree of arbitrariness in classifying earthquakes. When a 
relative big earthquake occurs, it can be either a foreshock 
of a larger subsequent event, an aftershock of a preceding 
large quake or the mainshock itself. Thus there is no way 
to differentiate these events in real time and many algo- 
rithms and methods are useless from a practical point of 
view as they perform only a posteriori. However, regard- 
less of how a big event may be classified, the forecasting 
of large quakes is the main goal to be achieved. This is the 
type of challenge that we face here. 

In this paper, we analyze a dissipative fibre bundle 
model of fracture with only one parameter that character- 
izes the heterogeneity degree of the system J2D3 - Notice- 
ably, the model exhibits two very distinct failure regimes 
obeying different avalanche size distributions. This allows 
the study of the occurrence of large avalanches, for the two 
phases, within the same model. In one of these regimes, the 
avalanche size distribution shows a power law regime for 
small events separated by a gap from the larger events, 
which are of the order of the system size. This regime 
is reminiscent of the characteristic earthquake behavior 
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found in some seismic faults [HIES- On the other hand, 
in the second regime, the system self-organizes into a crit- 
ical state characterized by a unique power law avalanche 
distribution similar to the GR law. Following the Molchan 
method [231 we compute the error diagrams that quantify 
the forecasting of large cascade events and analyze their 
correlations in time. Our results could be potentially ap- 
plied to the study of earthquake occurrence. 



2 Failure Model 

The basic ingredients and main features of the model used 
here are as follow |2(J| : 

1. A set of N elements is located on the sites of a sup- 
porting square lattice of linear length L. Each element 
represents a fibre, or in general terms, a small volume 
fraction of the material under study. 

2. To each element i, one assigns a random threshold 
strength ai th taken from a probability distribution. 
The threshold values represent the maximum load each 
of the elements can support before failure occurs. We 
use henceforth a Weibull distribution P(uth) = 1 — 
e -(f t f l /o'o) p to assign the failure thresholds assuming a 
reference load oq = 1. p is called the Weibull index; 
the bigger it is, the narrower the range of threshold 
values is |20"|. 

3. At each time step, the load on the system is increased 
by calculating the minimum load required for one ele- 
ment to break and adding this amount to all the ele- 
ments, i.e., <Tj — > <7i + where v — min{ai th — a{\. 

4. Once the load acting on an element surpasses its failure 
threshold, it fails and the load it bears is equally redis- 
tributed among all the elements of the set irrespective 
of their state. This amounts to a global fibre bundle 
version where long range interactions are assumed |24| . 

5. The fraction of load that would correspond to already 
failed elements is simply dissipated so that the model 
is nonconservative. Moreover, the degree of dissipation 
depends on both the avalanche size and its internal 
structure |20| . 

6. The redistribution of load to non-broken fibres may 
trigger more failure events. In such a case, step 4 is 
repeated until the system reaches a state where all the 
elements support a load below their respective thresh- 
old values. 

7. When an avalanche has come to an end, the external 
driving, step 3, is repeated and the elements that have 
been broken in the preceding avalanche are healed as- 
signing to them new random strengths and a load equal 
to zero. 

The model is thus constructed by implicitly taking 
into account some general properties of any fracture pro- 
cess. Two ingredients of our model are worth mention- 
ing. First, we assume dissipation through broken elements. 
This amounts to consider that the energy dissipated dur- 
ing a given fracture event is proportional to the size and 
the duration of that event. Besides, we incorporate healing 
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Fig. 1. Distributions of avalanche sizes for the two failure 
regimes that characterize the model. The upper part corre- 
sponds to a system that can be regarded as fairly homogeneous 
(p = 4) while the lower figure is obtained for a more heteroge- 
neous system made up of N = 50000 fibres (p = 2). The solid 
lines have slopes —2.5 and —1.5 respectively. After |20j. 



only after an avalanche has ended because when the ongo- 
ing avalanche is developing, the elements have no time for 
healing and broken regions should not accumulate stress 
in that time interval. We additionally note that the results 
presented here are valid, strictly speaking, in the limiting 
case of very slow driving. 

There are two different time scales in the model 20 . 
One refers to the external driving of the system and defines 
the natural time interval between failure events. We mea- 
sure the time elapsed until a given event k as the sum of 
the incremental amounts of load added at every external 
loading of the system prior to the event k. This defini- 
tion of time assumes that the rate at which the system is 
loaded can be considered constant, and hence that v is lin- 
early proportional to the time elapsed since the last event. 
This approximation holds, for example, in seismology. The 
second time scale of the process quantifies the avalanche 
lifetimes, and it is measured as lattice updates. That is, 
we consider the avalanches to be instantaneous with re- 
spect to the external driving and separated by an interval 
of time whose length is proportional to the amount of en- 
ergy added to the system during the last external input. 

The iteration of the rules sketched above leads the sys- 
tem to two different behaviors depending on the level of 
heterogeneity assumed when assigning the threshold val- 
ues j20] • Figure ^ shows the distributions of avalanche 
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Fig. 2. Error diagrams for the two failure regimes. Open cir- 
cles represent the trade-off between f a and f e when the system 
can be regarded as homogeneous (p = 4 and N = 100 ele- 
ments). Full circles have been obtained for a system made up 
of iV = 10 4 fibres and p = 2; the dashed line shows the effects of 
reducing the target threshold, s c ; the solid line corresponds to 
the case in which avalanches of any size are forecasted (equiv- 
alent to a random process). 

sizes, .P(s), for the two failure regimes of the model. The 
upper panel shows the pattern typical of characteristic 
quakes which are distributed quasi-periodically in time. 
This regime, where there are large fluctuations of the load 
on the system, is obtained when the bundle of fibres is ho- 
mogeneous (large values of the Weibull index p), that is, 
when the deviation of the thresholds from a mean value is 
not large. Note that as the system size is increased, large 
events are better resolved and there are no intermediate 
avalanches. 

A richer and more complex behavior is found when 
the system departs from homogeneity and the bundle of 
fibres gets heterogeneous (small values of p). In this case, 
irrespective of the system size, a power law distribution 
for the avalanche sizes is obtained. This is depicted in 
the lower panel of Figure ^ for p = 2 and N — 50000 fi- 
bres. As we shall see in the next section, the quasi-periodic 
nature of the first regime makes the forecasting of large 
avalanches feasible and eventually, for large system sizes 
(strictly speaking, in the thermodynamic limit), one can 
always forecast such events with no errors. This is not any- 
more the case when the system can be regarded as hetero- 
geneous, as non-trivial temporal correlations and cluster- 
ing of large events come into play making the forecasting 
of catastrophic avalanches a hard task. 

3 Error Diagrams and Temporal Correlations 

In order to inspect whether or not we would be able to pre- 
dict large events, we first compute the errors diagrams for 
both failure regimes. They are a quantitative measure of 
the success in forecasting the occurrence of a given event 
and were introduced in seismology several years ago by 
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Fig. 3. Time sequence of events, in a model with p = 2, with 
sizes equal to or greater than s. The time intervals have been 
rescaled so that the density of points is the same in both series. 
The system consist of N = 50000 elements. Note the clustering 
of large quakes (s > 10 3 ). 



Molchan 23 with the aim of rigorously evaluating differ- 
ent algorithms for large earthquakes prediction, and sub- 
sequently used by other authors 19,27, 25ll29| . 

An error diagram requires the identification of a tar- 
get to be forecasted using alarms. Then, one plots the 
trade-off between the fraction of target events that were 
not predicted, f e , and the fraction of time that the alarm 
was active, /„. Ideally, one would like to get the smallest 
number of failures to predict at the minimum duration of 
the alarms. Points close to the origin of an error diagram 
represent the best predictions. Hence, the efficiency of a 
given forecasting algorithm is measured by how fast f e 
decreases when increasing f a . In practice, one switches on 
the alarm at a given time interval r after the preceding 
event. If r = 0, the fraction of failures to predict is zero, 
but the fraction of alarm time is 1. On the contrary, if 
r is equal to (or greater than) the time interval between 
two successive events, f a = but f e = 1. In this way, the 
result of a completely random (Poissonian) process is a 
diagonal line. 

Let us now construct the error diagrams for each of the 
regimes of the fracture model under study. First, we de- 
fine the targets as the large avalanches whose sizes exceed 
a given threshold s c and fix the time interval r at which 
the alarms are switched on after every large failure event. 
By varying the length of r and computing the fractions 
of failures to forecast and of alarm times, f e (r) and / a (r) 
respectively, the error diagrams shown in Fig. are ob- 
tained, where the dependency on r has been eliminated. 
The figure clearly illustrates how different with respect to 
successful forecast the two failure regimes are. Open cir- 
cles represent the results of the algorithm when applied 
to the prediction of characteristic events in the parame- 
ter region when the system is homogeneous. Here, target 
events occur quasi-periodically and thus setting r close to 
the mean time interval between large avalanches ensures a 
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inter-occurrence time (At) 



Fig. 4. Distributions of inter-occurrence time At for a system 
of N = 10 4 elements and p = 2. The upper figure is in linear- 
log scale so that the distribution of At follows an exponential 
law when s c is small (s c > 10). The lower panel corresponds 
to s c equal to or greater than 10 2 so that only large avalanches 
are taken into account. In this case, the inter-occurrence times 
are power law distributed. 



successful prediction with a minimum of alarm time. For 
example, one can choose r so that f a w 0.3 and f e « 0.2, 
a rather good forecasting. Besides, increasing the system 
size improves the prediction of characteristic events since 
there appears a gap between small and large events mak- 
ing the definition of targets less noisy. In this case, the 
error diagram approaches a delta function at the origin as 
large failure events occur almost periodically in time. 

On the contrary, when the system gets heterogeneous, 
in the error diagram f a is always above f e (full circles). 
Here, we consider an avalanche as a target when its size 
falls on the edge of the distribution P(s), where finite 
size effects appear. Reducing the target thresholds move 
the error diagram close to the diagonal line. In fact, the 
occurrence of an event of any size is a random process as 
can be observed in the figure when the threshold is set 
equal to 1. 



1.0 




Fig. 5. Error diagrams of the BTW model for different values 
of the target thresholds. Note that irrespective of s c , the fore- 
cast is completely random in agreement with the fact that the 
BTW model does not show any kind of temporal correlations. 
The system size is N — 75 x 75. 

Another interesting information that can be extracted 
from the error diagrams for the second regime is the exis- 
tence of some degree of non-trivial temporal correlations. 
The departure from a random process when the target 
threshold is increased indicates that large events show 
clustering in time. This is clearly appreciated in Fig. |21 
where two sequences of quakes for a system made up of 
A" = 50000 elements and p = 2 has been represented. The 
upper series is the occurrence of events with sizes equal to 
or greater than 10 2 while the lower sequence is obtained 
when s c is equal to 10 3 (upper limit for events not affected 
by finite-size effects). 

We also present in Fig. 0] the distribution of inter- 
occurrence time At for two values of the target thresh- 
olds. We see that when s c is increased so that only large 
avalanches are taken into account, the distribution D {At) 
radically changes from an exponential decay (typical of 
a random process) to a decaying power law. This is in 
agreement with the error diagrams previously shown and 
provides further evidence that small avalanches are not 
correlated while large ones are. Without providing a fig- 
ure, we point out that the distribution of inter-occurrence 
time for the first failure regime is a well-behaved function 
with a clear mean value and a dispersion that decreases 
as the system gets larger. 

4 Discussion and Conclusions 

The results shown in the preceding section indicate that 
the clustering of large avalanches makes the prediction 
of such events a hard task. The shape of the error di- 
agrams and the high rate of prediction failures can be 
understood by noticing that the average inter-occurrence 
time between all large events is not a useful reference in 
this case. Two other quantities seem to be more funda- 
mental: the inter-cluster average time defined as the time 
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interval between clusters of large avalanches and the mean 
intra- cluster time defined as the time between large events 
within the same cluster. Since the former is much larger 
that the latter, the fraction of failures to predict is raised 
when r is large. On the other hand, setting r too small to 
match the intra- cluster time produces a significant incre- 
ment of f a - Thus both phenomena contribute to the error 
diagrams in such a way that the forecasting is worse than 
a random prediction. 

The clustering of large events is due to the way in 
which energy is dissipated and not to the power-law na- 
ture of the avalanche size distribution. The average inter- 
occurrence time between large quakes can be considered 
as a measure of the buildup time for correlations in the 
system. When a series of big cascades occurs, the system 
ends up in a state where most of the elements are unloaded 
and a significant input of energy (or a large time interval) 
is needed in order to start another series. In other words, 
the way in which energy is dissipated creates stress corre- 
lations giving rise on its turn to the observed clustering. 
This does not hold anymore when avalanches of smaller 
sizes are considered. In this case, the occurrence of a small 
avalanche does not significantly alter the state of the sys- 
tem. A small event can not influence the occurrence of any 
event at short times. 

The same kind of temporal correlations is seen in other 
models displaying self-organized criticality. This is the case 
of the Olami, Feder and Christensen (OFC) spring-block 
model for earthquakes [7||S]. They found some years ago 
that non-trivial temporal correlations of the sort obtained 
here show up when the model is made nonconservative. 
By calculating the coefficient of variation, it was shown 
that large earthquakes show clustering while small ones 
are uncorrelated in time. Besides, when the model is con- 
servative, neither spatial nor temporal correlations were 
found 8 . The fact that dissipation appears to be a basic 
ingredient for a model to show correlations in time and/or 
space is corroborated in Fig. |SJ where the error diagrams 
for the original Bak, Tang and Wiesenfel (BTW) model 
[131114) are depicted. As can be clearly seen in the figure, 
the model can not be distinguished from a random pro- 
cess with regard to its predictability. Changing the tar- 
get thresholds does not improve the forecasting. Thus, in 
this model, small avalanches and large ones are completely 
equivalent and all correlations are ruled out. 

We would like to remark that the results obtained in 
this paper may be applied to the study of earthquake 
occurrence. One of the few well-known facts in seismol- 
ogy is the clustering of earthquakes in time and space. 
While small earthquakes seem to be uncorrelated |3Uj . 
large earthquakes display strong clustering |31ll32j . On the 
other hand, by tuning a single parameter in our model, 
one moves from a characteristic earthquake scenario |33| 
to a regime where the avalanche sizes are power law dis- 
tributed (GR law) . The possibility of exploring such differ- 
ent regimes within the same model is another motivation 
for further research in this direction |34j . 

Summing up, we have analyzed a fracture model with 
non-trivial temporal correlations. The results presented 



seem to support the idea that models aimed at modeling 
real earthquakes should incorporate some kind of dissipa- 
tion and self-organization. On the other hand, the clus- 
tering of large quakes does not help to improve the pre- 
dictability of devastating events since two relevant time 
scales come into play. It might be interesting to compute 
the error diagrams with a tunable alarm time in such a 
way that one can better resolve series of large events and 
periods of stasis. 

Y. M. acknowledges financial support from the Secretaria de 
Estado de Education y Universidades (Spain, SB2000-0357). 
M. V-P. is supported by the Ph.D grant B037/2001 of the 
DGA. This work has been partially supported by the Spanish 
DGICYT project BFM2002-01798. 



References 

1. H. J. Herrmann, and S. Roux, (Eds), Statistical Models for 
the Fracture of Disordered Media (North Holland, Amster- 
dam, 1990), and references therein. 

2. D. L. Turcotte, Fractals and Chaos in Geology and Geo- 
physics. 2nd ed (Cambridge University Press, Cambridge, 
New York, 1997). 

3. D. Sornette, Critical Phenomena in Natural Sciences, 
(Springer- Verlag, New York, 2000) . 

4. L. de Arcangelis, S. Redner, and H. J. Herrmann, J. Phys. 
(Paris) 46, L585 (1985). 

5. S. Roux, and E. Guyon, J. Phys. (France) Lett. 46, L999 
(1985). 

6. R. Burridgeq, and L. Knopoff, Bull. Seismol. Soc. Am. 57, 
341 (1967). 

7. Z. Olami, H. J. Feder, and K. Christensen, Phys. Rev. Lett. 
68, 1244 (1992). 

8. K. Christensen, and Z. Olami, J. Geophys. Res. 97, 8729 
(1992). 

9. H. E. Daniels, Proc. R. Soc. London, Ser. A183,404 (1945). 

10. W. I. Newman, D. L. Turcotte, and A. M. Gabrielov, Phys. 
Rev. £52, 4827 (1995). 

11. P. C. Hemmer,and A. Hansen, J. Appl. Mech. 59, 909 
(1992). 

12. Y. Moreno, J. B. Gomez, and A. F. Pacheco, Phys. Rev. 
Lett. 85, 2865 (2000). 

13. P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 
381 (1987). 

14. P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A. 38, 364 

(1988) . 

15. G. Gutenberg, and C. F. Richter, Ann. Geophys. 9, 1 
(1956). 

16. P. Bak, and C. Tang, J. Geophys. Res. 94, 15635 (1989). 

17. D. Sornette, and A. Sornette, Europhys. Lett. 9, 197 

(1989) . 

18. K. Ito, and M. Matsuzaki, J. Geophys. Res. 95, 6853 

(1990) . 

19. V. Keilis-Borok, Annu. Rev. Earth Planet. Sci. 30, 1 
(2002). 

20. Y. Moreno, J. B. Gomez, and A. F. Pacheco, Physica A 
274, 400 (1999). 

21. S. G. Wesnousky, Bull Seismol. Soc. Am. 84, 1940 (1994). 

22. K. Sieh, Proc. Natl. Acad. Sci. USA 93, 3764 (1996). 



6 



Yamir Moreno et al.: Error Diagrams and Temporal Correlations. 



23. G. M. Molchan, Pure Appl. Geophys. 149, 233 (1997). 

24. This choice may seem arbitrary. While the equal load shar- 
ing is the simplest load transfer rule one may consider, 
it also gives more accurate results when modeling elastic 
forces than other load transfer schemes, at least within the 
FBM context. However, one can in principle consider gen- 
eral load transfer schemes as the one introduced in |25|. 
There, it is shown that the behavior of the system can be 
regarded as global provided that the interactions follows 
a power law of the form r~"' with 7 < 2. On the other 
hand, similar results have been obtained when using local 
models |26|. 

25. R. C. Hidalgo, Y. Moreno, F. Kun, and H. J. Herrmann, 
Phys. Rev. E 65, 046148 (2002). 

26. Z. Olami, H. J. Feder, and K. Christensen, Phys. Rev. Lett. 
68, 1244 (1992). 

27. W. I. Newman, and D. L. Turcotte, Nonhn. Proc. Geophys. 
9, 1 (2002). 

28. A. Helmstetter, and D. Sornette, e-print 
cond- mat/0208597 (2002). 

29. M. Vazquez-Prada, A. Gonzalez, J. B. Gomez, and A. F. 
Pacheco, unpublished (2003). 

30. A. C. Johnston, and S. J. Nava, J. Geophys. Res. 90, 6737 
(1985). 

31. Y. Y. Kagan, and D. D. Jackson, Geophys. J. Int. 104, 
117 (1991). 

32. M. S. Mega, P. Allegrini, P. Grigolini, V. Latora, L. 
Palatella, A. Rapisarda, and S. Vinciguerra, Phys. Rev. 
Lett, in press. Also e-print cond- mat/0208597 (2002). 

33. M. Vazquez-Prada, A. Gonzalez, J. B. Gomez, and A. F. 
Pacheco, Nonlin. Proc. Geophys. 9, 513 (2002). 

34. Models able to pass from the GR to a characteristic earth- 
quake behavior are not found very often in the literature. 
One of such models was proposed in |35j . where the con- 
cept of configurational entropy is used to argue the ex- 
istence of the two regimes. Our model also suggests that 
these two behaviors can be recognized when the degree of 
heterogeneity of the faults under study are different. 

35. K. Dahmen, D. Ertas, and Y. Ben-Zion, Phys. Rev. E 58, 
1494 (1998). 



